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Bottomonium S-wave states were studied using lattice NRQCD. Masses of ground and excited 
states were calculated using multiexponential fitting to a set of correlation functions constructed 
using both local and wavefunction-smeared operators. Three-point functions for Ml transitions 
between vector and pseudoscalar states were computed. Robust signals for transitions involving the 
first two excited states were obtained. The qualitative features of the transition matrix elements 
are in agreement with expectations. The calculated values of matrix elements for T(25) and T(35) 
decay are considerably larger than values inferred from measured decay widths. 



I. INTRODUCTION 

The bottomonium T was discovered in 1977 []], [2j 
and, remarkably, it took 30 years before its pseudoscalar 
partner rji, was observed^, |j|. The measurement of the 
branching ratio for radiative decay of T(2S) and T(35) 
to r\b presents an opportunity to test calculational meth- 
ods where the decay amplitude depends entirely on small 
effects: spin-dependent interactions, recoil and relativis- 
tic corrections. In this paper we present a first pass at 
the calculation of the excited Upsilon radiative decays 
using lattice NRQCD@. 

To understand the challenge of excited Upsilon decays 
it is useful to consider first the amplitude for the magnetic 
dipole (Ml) transition between vector and pseudoscalar 
states in the nonrelativistic quark mo 
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where R n ( r ) is the radial wavefunction for the S-wave 
state with principal quantum number n and q is the pho- 
ton momentum. In a transition between states with the 
same principal quantum numbers n' — n, for example 
ground state to ground state, the radial wavefunctions 
of the vector and pseudoscalar mesons are very similar. 
The overlap integral is close to 1. However, when n' =/= n 
the wavefunctions are orthogonal in the extreme nonrel- 
ativistic limit. For these so-called hindered transitions 
the amplitude is highly suppressed and depends on the 
interplay of small effects coming from spin-dependent in- 
teractions, effects of recoil and relativistic corrections [6- 

The use of lattice QCD methods to calculate the am- 
plitude for vector meson radiative decays was suggested 
long ago [j| [nj. Recently, interest in this application 
of lattice QCD has been revived and charmonium has 
been studied in detail fill— Tl~3l| . A number of different 
ground-state to ground-state transition amplitudes have 
been calculated involving not just S-wave but also P-wave 
states [ll|, [HI . In Ref. [l2j radiative decays of excited 
charmonium states were also considered. Excited states 



appear as nonleading contributions to lattice QCD meson 
correlation functions. This, combined with the suppres- 
sion of the hindered Ml amplitude, makes it a challenge 
to achieve good statistical accuracy for these decays (see 
Table III in Ref. 

The application of lattice QCD to excited states is now 
a very active research area. A primary goal of this work 
is to see how well we can extract the excited state sig- 
nal buried under the dominant ground-state contribu- 
tion. One way to deal with this problem is to use a vari- 
ational method [id] . [l5| with an appropriate (and large) 
set of basis operators. An alternative which has been 
applied successfully to the calculation of the spectrum 
of bottomonium is to use constrained multiexponential 
lilting. l(i| to get the subdominant contributions. This 
method can work well if the lattice simulation data have 
high statistical precision (see, for example, Ref. |l7|). 
The gauge field ensemble used in this study is quite small, 
only 192 configurations, but we reduce the statistical fluc- 
tuations by using multiple time sources per configuration 
and by employing a spatial wall source. This allows the 
extraction of robust signals for the 2S and 3S excited 
states. 

Section |TT] outlines the lattice QCD simulation. The 
gauge field configurations come from a 2+1 flavor dy- 
namical simulation and were provided by the PACS-CS 
Collaboration[l8j. The b quarks are described using a 
standard 0{v A ) lattice NRQCD action [!, [H, [20| with 
Landau link tadpole improvement. Two-point correla- 
tion functions of pseudoscalar and vector operators are 
discussed in Sec. Mil Two-point function fit parame- 
ters, simulation energies and overlap coefficients were 
obtained, and these are used unchanged in subsequent 
three-point function fits. As a check of the calculation the 
simulation energies for the lowest three states from our 
multiexponential fits were compared to the variational 
analysis that could be done with our limited basis set. 
Good consistency was obtained. As well, the overlap co- 
efficients for the lowest three Upsilon states were used 
to estimate leptonic decay widths with reasonable agree- 
ment with experiment. Three-point functions and tran- 
sition matrix element results are discussed in Sec. IIV1 A 
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summary is given in Sec. [V] 



II. LATTICE SIMULATION 

The lattice gauge field configurations were generated 
and made available by the PACS-CS Collaboration|18j . 
The 2+1 flavor dynamical simulation used the Iwasaki 
gauge field action (/3=1.90) and the clover- Wilson 
fermion action. We use only the ensemble which is 
nearest the light-quark physical point with pion mass 
156MeV. The light and strange hopping parameters are 
K u/d = 0.13781 and k s — 0.13640. The number of 
lattice points is 32 3 x64 and the lattice spacing a = 
0.0907(14)fm was determined by PACS-CSQ (along 
with the light and strange hopping parameters) using the 
pion. kaon and Q~ baryon masses as input, i.e., heavy- 
hadron input was not used in setting the scale. The num- 
ber of gauge field configurations used was 192 (out of 
198 available). For the tadpole-improved NRQCD calcu- 
lation the average link in Landau gauge was used. The 
numerical value was estimated to be 0.8463. 

The heavy quark is described using lattice NRQCD H, 
IT9L [20| . The exact form of the Hamiltonian may be found 
in the Appendix of Ref. [21]. Terms up to 0(v 4 ) are kept 
in the nonrelativistic expansion which in the notation of 
[2l| means a = 1 for i < 6 and c, = for i > 7. The 
b-quark bare mass was determined by fitting the kinetic 
mass to the observed mass of rjt and takes the value 1.945 
in lattice units. The stability parameter n appearing in 
the Hamiltonian was taken to be 4 in line with Ref. [17j . 

The simplest operators to use to describe the pseu- 
doscalar and vector states are the local ones which take 
the form 0(x) = x(x)Tip(x) where ip(x) and x( x ) are 
nonrelativistic quark and antiquark fields with T equal 
1 (a) for pseudoscalar (vector) mesons. To calculate 
ground-state properties a smearing of the local operators, 
such as Jacobi smearing [2^, is often used to damp out the 
high-energy modes created by local operators. However, 
for investigating excited states it is more advantageous to 
include operators which suppress the ground state. For 
this, wavefunction smearing[19] is useful. The smeared 
operator takes the form 0{x) = ^2 y x(x)T(f>(x — y)tp(y) 
where an effective smearing function was found to be(l7| 



(r) = (1 - r/(2 a o)K r/(2ao) 



(2) 



which has the profile of the Coulomb S-wave first-excited 
state wavefunction. The parameter ao was taken to be 
1.4(lattice units). In addition to the smeared operator a 
doubly-smeared operator where the wavefunction smear- 
ing was applied to both quark and antiquark fields was 
included. Our complete set of meson operators consisted 
of three types: local (1), smeared (s) and doubly-smeared 
(d). In order to use the nonlocal smeared operators with- 
out gauge links connecting the quark and antiquark the 
gauge field configurations were fixed to Coulomb gauge. 
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Figure 1: Kinetic energy of the pseudoscalar meson at b-quark 
bare mass Ma = 1.945 versus momentum squared. The line 
is a fit with Eq. ©. The x 2 /d.o.f. is 0.2. 



III. TWO-POINT FUNCTIONS 

In NRQCD the correlators of the meson operators do 
not give the hadron mass directly. The simulation energy 
extracted from the zero-momentum correlation function 
must be combined with the renormalized quark mass and 
energy shift to get the meson mass. Alternatively the ki- 
netic energy can be used to determine the meson mass. 
This is the method used here for tuning the b-quark mass 
to reproduce the mass of rji,. Correlators of the local pseu- 
doscalar meson operator projected onto different values 
of momentum were calculated. Using the relation 



E(p) - £7(0) = (p 2 + M 2 ) 



1/2 



M 



(3) 



the hadron mass Mo can determined. Using this method 
we arrive at a bare quark mass value of 1.945(4). The 
quoted error in this value reflects the uncertainty in the 
fit determining the kinetic mass. As well, there is a 1.5% 
uncertainty due to the uncertainty of the lattice spac- 
ing determination. The kinetic energy and the fit that 
determines Mo at our nominal bare mass are shown in 
Fig. [TJ The momenta used were (0,0,0), (1,0,0), (1,1,0) 
and (1,1,1) in units of 2n/La where La is the spatial 
extent of that lattice. 

Using the determined value of the bare b-quark mass 
two-point correlation functions (including cross correla- 
tors) of the three operator types l,s,d were calculated for 
pseudoscalar and vector channels. Our lattice has 64 time 
sites but it is not useful to construct correlators over the 
entire time extent. The correlators are limited to max- 
imum time separation t — t s of 27. Since the maximum 
time separation is considerably smaller than the lattice 
time extent and the nonrelativistic propagators depend 
only on the gauge field links on time slices between source 
and sink it is very effective to use multiple time sources 



3 



o 10 
c 



lio- 3 

o 
U 

io- 4 



r^^ti 1 i 1 i 1 i 1 r 









Sis 




Sid 




Figure 2: Zero-momentum vector correlation functions for dif- 
ferent operator combinations. Symbols are simulation values 
and lines are the result of a fit with ten exponential terms. 
Except for some points with the ss operator combination, sta- 
tistical errors are smaller than the symbols. 



on each gauge configuration. Sixteen sources, uniformly 
distributed in time, were used in this calculation. To fur- 
ther reduce statistical fluctuations a random wall source 
(see, for example, [23]) was used. In addition to zero- 
momentum pseudoscalar and vector meson correlators, 
pseudoscalar correlators with momentum corresponding 
to (1,0,0) and (2,0,0) were calculated. These are needed 
for the analysis of the three-point functions. 

Using subscripts o and d (o, d = {I, s, d}) to de- 
note source and sink operators, the correlation functions 
goo'(t) = (0 '(t)Ol(t s )) (fixed source time t s ) are fit with 
N time-dependent exponentials 



9oo'{t) = ^ c -(n)c (n)e 



-E n (t-t s 



(4) 



The constrained multiexponential fitting method[lg, [17j 
was used. All time points (except the source) were in- 
cluded. Using only loose constraints fits are very stable 
even with a large number of exponential terms. Figure 
[2] shows a representative sample of correlation functions 
and fits, in this case for the zero-momentum vector chan- 
nel with a simultaneous fit using 10 terms. 

The lowest four simulation energies for zero- 
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Figure 3: The simulation energies of the four lowest pseu- 
doscalar and vector states from a multiexponential fit to the 
two-point functions with 6, 8, and 10 terms. 
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Figure 4: Mass difference between the Upsilon ground state 
and the first two excited states using results of 6, 8 and 10 
term multiexponential fits. The dashed lines show experimen- 
tal values using data from [24]. 



momentum pseudoscalar and vector channels are shown 
in Fig. [3] for fits with 6, 8 and 10 terms. The lowest 
three states, which are of interest for our three-point 
function calculation, are quite robust. Differences of the 
zero-momentum simulation energies are just mass differ- 
ences and these are shown for the lowest Upsilon states 
in Fig. 2J The results are in reasonable accord with ex- 
perimental values jiH. 

The mass difference between T and rjb was also cal- 
culated. For the ground states, our result is 56(l)MeV 
where the error is dominated by the uncertainty in the 
lattice spacing. This value is essentially independent of 
which fit is used and consistent with values obtained 
by others [20|, EH using the 0{v 4 ) lattice NRQCD ac- 
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Figure 5: The effective mass of the eigenvalues from a vari- 
ational analysis of the correlator matrix for the pseudoscalar 
meson. The lines show the simulation energies for the lowest 
three states from the 10-term multiexponential fit. 



Figure 6: The effective mass of the eigenvalues from a varia- 
tional analysis of the correlator matrix for the vector meson. 
The lines show the simulation energies for the lowest three 
states from the 10-term multiexponential fit. 



tion. It is somewhat smaller than the PDG average 
69.8±2.8MeV of the experimentally observed values[3, 
|26| . It is a common feature for lattice simulations of 
heavy quark systems to underestimate the spin splitting 
and some issues have been discussed in the context of 
lattice NRQCD in recent studies [U HI- 

A popular way to determine excited state energies is 
the variational method (Til [l5| (for an extensive recent 
discussion see [29]). The correlator matrix is diagonal- 
ized at each time and the time-dependent eigenvalues 
then give an optimal estimate of the time evolution of in- 
dividual states. The evolution is calculated with respect 
to some reference time to > t s which should be chosen 
large enough so that the number of basis operators is 
comparable to the number of contributing states. Our 
operator set is too small to use this method effectively 
but it is of interest nonetheless to compare this method 
with the results of the multiexponential fit. Figs. [5] and 
[fj] show the effective mass plots for the eigenvalues A*, (t) 
which are solutions of the generalized eigenvalue problem 
(with to = 4) 
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(5) 



where fk(t) is the eigenvector. For higher states only 
a limited number of time steps are available before the 
effective mass degenerates into noise. The lines on the 
plots show the simulation energies from the 10-term mul- 
tiexponential fit. For the ground and first-excited states, 
where a meaningful comparison is possible, there is com- 
plete consistency. 

The overlap coefficients in the vector channel provide 
another test of the calculation. They can be used to de- 
termine the partial width for Upsilon states to decay into 



Table I: Decay amplitude and partial width for Upsilon lep- 
tonic decay. The exp erimental values T e x P are from the Par- 
ticle Data Group[24]. 



State 



,3/2 



tt(0) TjkeV] r exp [keV] 

T(1S) 0.18418(7) 1.16(2) 1.34(2) 

T(2S) 0.1397(33) 0.595(28) 0.612(11) 

T(3S) 0.156(24) 0.70(21) 0.443(8) 



lepton pairs and this can be compared to experimental 
values. The decay width can be expressed m terms of the 
wavefunction at the origin V I'„(0) as (see [20]) 



T(T(nS) 
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where a is the electromagnetic coupling constant and 
the matching factor Z matc h relates the lattice vector cur- 
rent to the renormalized continuum current. With non- 
relativistic normalization of states, ^(0) is related to 
the overlap coefficient of the local vector operator by 
\l/(0) = q/ The results from our calculation are shown 
in Table HI The matching factor Z matc h has not been cal- 
culated for the version of lattice NRQCD that we use so 
the leading order value of 1 has been assumed. Hart et 
al. [30] have computed that matching coefficient for lattice 
NRQCD with stability parameter equal 2. Using their 
result as a guide we might expect effects from match- 
ing of about 5% but without an actual calculation it is 
not possible to say with certainty which way they would 
go. Given this state of the calculation, consistency with 
experiment is reasonable. 
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IV. THREE-POINT FUNCTIONS 

Three-point functions describing the vector to pseu- 
doscalar transition induced by a current operator inser- 
tion are constructed using a sequential source method. 
The transition operator is taken here to be just the lead- 
ing nonrelativistic operator a which acting on a quark (or 
antiquark) converts a vector state to a pseudoscalar (or 
vice versa). Starting with a vector (pseudoscalar) source 
at t s the quark propagator is evolved to some maximum 
time T at which a pseudoscalar (vector) operator is ap- 
plied. This quantity is then evolved backward in time. At 
intermediate times t s < t' < T the current operator is in- 
serted and evolution is continued to complete the quark 
antiquark loop at the source. Appropriate momentum 
projections are applied at the source, sink and current 
insertion to ensure momentum conservation. The vec- 
tor meson is always projected to have zero momentum; 
the pseudoscalar recoils against the momentum carried 
by the current. 

With a vector operator at the source and pseudoscalar 
at the sink the three-point function is expected to have 
the form 



Table II: Three-point matrix elements from simultaneous fits 
to N c f correlation functions and with different numbers of 
parameters at zero recoil momentum. 



xe -E { ^(t'-t s ) e -E^(T~t') 



(7) 



where the subscripts o, d indicate the type of smearing 
used (l,s or d). The overlap coefficients and simulation 
energies are the same ones that appear in the two-point 
function but now have a superscript attached to distin- 
guish between vector and pseudoscalar states. The quan- 
tity A^, is the matrix element of the transition op- 
erator between the vector state n and the pseudoscalar 
state n' . This is identified with the wavefunction overlap 
appearing in ([1]). The three-point function with pseu- 
doscalar source and vector sink has the same form with 
V and P labels reversed. The matrix elements Al^ 

are related to those appearing in (0 by A!j^P = ^Xf ■ 
The matrix elements can be determined by fitting the t' 
dependence of the three-point function for a fixed T. 

If the spin-dependent interaction terms in the NRQCD 
Hamiltonian, which are nonleading in the nonrelativistic 
expansion, were omitted and pseudoscalar meson recoil 
momentum set to zero, the three-point functions would 
be independent of t' and numerically equal to the two- 
point functions at time separation T — t s for all T . The 
solution for the matrix elements would be trivial A^S — 
8 nn i , the same as for the wavefunction overlap ([1} in the 
extreme nonrelativistic limit [7|. 

Three-point functions were computed for three val- 
ues of pseudoscalar recoil momentum corresponding to 
(0,0,0), (1,0,0) and (2,0,0) and for two source-sink time 
separations T — t s equal to 27 and 19. All combinations 
of operator types l,s,d were calculated but in the analy- 
sis only the combinations 11, Is, si, Id, dl, ss were used. 
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a(pv) 




a(VP) 


A (VP) 






T 


- t s = 19 






10 


0.916(2) 


-0.043(7) 


-0.069(6) 0.090(7) 


0.052(5) 




10 


0.915(2) 


-0.068(2) 


-0.050(4) 0.072(4) 


0.065(3) 


1.11(31) 


12 


0.915(2) 


-0.068(3) 


-0.050(4) 0.071(4) 


0.065(3) 


1.11(23) 






T 


- t s = 27 






10 


0.916(2) 


-0.062(7) 


-0.056(7) 0.075(7) 


0.059(6) 




10 


0.916(2) 


-0.068(3) 


-0.050(6) 0.071(3) 


0.062(4) 


2.1(2.2) 


12 


0.916(2) 


-0.068(3) 


-0.051(6) 0.071(4) 


0.062(4) 


1.9(1.8) 



Fits were done both including and excluding the ss op- 
erator combination. The other combinations were noisy 
and not useful in pulling out the excited to ground-state 
transitions that are of interest here. 

The two-point functions were fit with many exponen- 
tial terms in order to get a stable result for the lowest few 
states. For the three-point function with large source- 
sink separation the contribution of the high-lying states 
is highly suppressed and neglected in the fit of the t' de- 
pendence. Only the lowest three states were considered 
and the nn' combinations included in the sums in ([7]) 
were 11, 12, 21, 13, 31, 22. As a test of the robustness 
of the results some five-parameter fits, excluding the 22 
term, were done. The time fit range was taken to be 
t s ■ + 2 < t' < T — 2. The fits using two-point function 
parameters, overlap coefficients and simulation energies, 
determined using ten terms are given here. Using eight- 
term two-point function parameters gives essentially the 
same values. The six-term two-point function parameters 
lead to slightly different results but we do not consider six 
terms to be sufficient for the two-point function fit. The 
determination of the matrix elements is done by a simul- 
taneous fit to a set of three-point functions. Statistical 
errors are estimated using a bootstrap analysis. Some 
representative three-point correlation functions and fits 
are given in Figs. and [5] for T — t s equal to 19 and in 
Figs. M and [TO] for T - t s equal to 27. 

The results for the three-point matrix elements are 
given in Tables |n] - IIVI for different values of recoil mo- 
mentum. For T — t s = 19 the 2 to 2 transition is clearly 
necessary to get results that are consistent with the larger 
time separation. For the excited to ground-state transi- 
tions, that are of primary interest, there is very good 
agreement between results using the shorter and longer 
time separation. 

The results from the T — t s — 19 analysis with a six 
parameter fit are plotted in Fig. [11] to [T3] as a function 
of momentum. In Fig. Q2] the matrix elements inferred 
from the measured T(25) and T(3S*) to r/b partial widths 
[H, 0| are also shown at the physical momentum for these 
decays. The features of the results are easily understood. 
The matrix elements for the T(15) to r]b(lS) and T(2S) 
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Figure 7: Three-point correlation functions with vector source 
and pseudoscalar sink at time separation T — t s equals 19 
for different operator combinations. Symbols are simulation 
values and the lines are the result of a simultaneous fit. The 
t' fit range is 4 to 17. 



Figure 8: Three-point correlation functions with pseudoscalar 
source and vector sink at time separation T — t a equals 19 
for different operator combinations. Symbols are simulation 
values and the lines are the result of a simultaneous fit. The 
t' fit range is 4 to 17. 



Table III: Three-point matrix elements from simultaneous fits 
to N c f correlation functions and with different numbers of 
parameters at one unit of recoil momentum. 
N ' a ^p) TTpvT TTpvT ~A~CyP) ~a(vp) A i vp ) 

T-t s = 19 

10 0.908(1) -0.042(8) -0.060(8) 0.095(7) 0.057(5) 

10 0.907(1) -0.062(6) -0.047(7) 0.079(4) 0.068(5) 0.92(27) 

12 0.907(1) -0.062(6) -0.047(7) 0.079(5) 0.067(5) 0.95(21) 

T - t s = 27 

10 0.908(2) 0.057(8) -0.052(9) 0.082(6) 0.063(6) 

10 0.907(2) 0.061(5) -0.048(8) 0.079(4) 0.066(6) 1.6(1.9) 

12 0.907(2) 0.061(5) -0.048(8) 0.079(5) 0.066(6) 1.6(1.5) 



Table IV: Three-point matrix elements from simultaneous fits 
to N c f correlation functions and with different numbers of 
parameters at two units of recoil momentum. 
N ' ~A<yP) JiPV) aTpv) TTFp) T(vp) a (vp) 

iV e/ ^11 ±\21 1221 1121 ±\S1 A 22 

T — t s — 19 

10 0.878(1) -0.010(6) -0.055(6) 0.116(7) 0.066(6) 
10 0.877(1) -0.030(4) -0.041(6) 0.101(5) 0.078(6) 1.01(25) 
12 0.877(1) -0.031(4) -0.041(6) 0.102(5) 0.078(6) 1.02(20) 
T — t s — 27 

10 0.878(1) -0.026(6) -0.041(8) 0.104(6) 0.066(8) 

10 0.878(2) -0.031(4) -0.037(8) 0.101(5) 0.070(6) 1.9(1.8) 

12 0.878(2) -0.029(5) -0.039(7) 0.100(5) 0.068(6) 1.0(1.6) 



to r]b(2S) transitions are close to 1 since the wavefunc- 
tions of the states involved are very similar. The T(IS') 
to r]b(lS) matrix element decreases only very slowly with 
recoil momentum which reflects the small size of bot- 
tomonium. The excited state to ground-state transitions 
have matrix elements that are small in magnitude due to 
near orthogonality of wavefunctions. The relative nega- 
tive sign of T — » 775 and rjb — > T reflects the fact that 
the dominant spin-dependent quark antiquark interac- 
tion acts with different sign in pseudoscalar and vector 
states. The recoil effect contributes positively to all tran- 



sitions (see, for example, [7|) which explains the momen- 
tum dependence. 

The calculated matrix elements are large compared 
to the empirical values inferred from the measured par- 
tial widths. However, there are a variety of improve- 
ments that are needed before definitive conclusions can 
be drawn. The lattice vector current operator has to be 
matched to the renormalized continuum current as dis- 
cussed in connection to Upsilon leptonic decay [5(3, [30j . 
Relativistic corrections are not incorporated into the 
transition operator and these are likely to be important 
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Figure 9: Three-point correlation functions with vector source 
and pseudoscalar sink at time separation T — t s equals 27 
for different operator combinations. Symbols are simulation 
values and the lines are the result of a simultaneous fit. The 
t' fit range is 4 to 25. 
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Figure 10: Three-point correlation functions with pseu- 
doscalar source and vector sink at time separation T — t s 
equals 27 for different operator combinations. Symbols are 
simulation values and the lines are the result of a simultane- 
ous fit. The t' fit range is 4 to 25. 



for the hindered excited state decays [7j. As well, one 
might ask about 0(v 6 ) terms j2?| and radiative correc- 
tions (beyond tadpole improvement) to spin-dependent 
interactions in the Hamiltonian which have been shown 
to have a noticeable effect on the Y—rjb mass splitting (28|. 
Finally, there is the question of continuum extrapolation 
which this calculation, done at a single lattice spacing, 
can not address. 

There are other systematics that we can not deal with 
quantitatively but which it is reasonable to think are 
small. Bottomonium has only heavy valence quarks so 
extrapolation to the physical point for up and down 
quarks comes in only through the influence of sea quarks 
on the gauge field. Since the simulation is done very near 
the physical point, with quarks which give a pion mass 
of 156MeV, it would be surprising if a simulation at the 
physical point would be much different. Finite volume 
can also lead to a significant systematic effect in lattice 
simulations but is unlikely to be the case here. The size 
of bottomonium is much smaller than the spatial lattice 
size so finite volume effects arise indirectly through light 
quarks in the sea. We do not have any estimates of this 
effect. However, finite volume effects have been studied 
for heavy-light mesons For our lattice size they are 
very small and it is reasonable to expect them to be even 
smaller for bottomonium. 



V. SUMMARY 

Bottomonium S-wave states were studied using lattice 
NRQCD focusing on the low-lying excited states. It was 
found that using a set of operators, including smeared 
operators which suppress the ground-state contribution 
to the correlation functions, robust results for the lowest 
few states could be obtained. Constrained multiexpo- 
nential fittingpjl was used for the analysis of two-point 
correlation functions. As a check, an analysis based on 
the variational methodfhU. [l5| was carried out. Where a 
meaningful comparison could be made, the two analysis 
methods gave consistent results. 

Mass differences between the Upsilon ground state and 
the first two excited states are in reasonable agreement 
with experimental values. The mass difference between 
T and r\b is not well reproduced by the calculation. Issues 
such as continuum extrapolation, higher order nonrela- 
tivistic terms and radiative corrections to the NRQCD 
Hamiltonian have to be considered. 

A primary goal of this study was to see if the highly 
suppressed matrix elements of excited state transitions 
could be extracted. Three-point functions for transitions 
from vector to pseudoscalar states with the leading non- 
relativistic Ml operator were calculated. Using overlap 
coefficients and simulation energies obtained from fitting 
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Figure 11: The matrix elements for decay of T to r/b with the 
same principal quantum number as a function of momentum. 



Figure 13: The matrix elements for decay of an excited rjb to 
the T ground state as a function of momentum. 
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Figure 12: The matrix elements for decay of an excited T to 
the rjb ground state as a function of momentum. The square 
symbols show the matrix element values inferred from the 
measured decay widths 0, 0]. 



the two-point functions, a simultaneous fit was done to 
sets of three-point functions. Transition matrix elements 
with reasonably small statistical errors could be obtained 
for a number of excited state decays. The results were 
very stable with respect to choice of two-point function 
parameters and the number of three-point functions and 
matrix elements included in the fit. 

The qualitative features of the calculated matrix el- 
ements are as expected. The matrix elements for the 
T(liS) to rjb(lS) and T(25) to r]b(2S) transitions are close 
to one. For states identified with different principal quan- 
tum numbers the transitions are highly suppressed. The 
relative negative sign of T — > rj b and rjb — > T matrix el- 
ements can be understood by considering perturbatively 
the effect of spin-dependent interactions. The qualitative 
momentum dependence is in accord with, for example, 
pNRQCDQ. 

Quantitatively the values obtained here for excited 
T to ground-state rjb matrix elements are considerably 
larger than the values inferred from the experimentally 
determined decay widths. These decays are dependent 
on the interplay of small effects and are likely to be 
sensitive to relativistic corrections to the transition op- 
erator. As well, the issues that enter into the T - rjb 
spin splittin g, e. g., relativistic corrections to the NRQCD 
Hamilt onian 12 7ll . operator matching [28j and continuum 
extrapolation[20] have to be dealt with to get definitive 
results. 
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